Clustering

load(paste0("./GenesCO2_", specie, ".RData"))


load("./normalized.count_At.RData")

# quantification file
data <- read.csv("quantifFiles/QuantifGenes.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, grepl("R", colnames(data))]

keep <- rowSums(data) >= 10
data <- data[keep, ]
group <- sapply(colnames(data), getLabel, with.rep = F)
colnames(data) <- sapply(colnames(data), getLabel)

specie = "At"
clusteredGenes <- clustering(sharedBy3, round(normalized.count, 0))

****************************************
coseq analysis: Poisson approach & none transformation
K = 2 to 12 
Use set.seed() prior to running coseq for reproducible results.
****************************************
Running g = 2 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0"
Running g = 3 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0"
Running g = 4 ...
[1] "Initialization: 1"
[1] "Log-like diff: 1.72263980857679e-10"
Running g = 5 ...
[1] "Initialization: 1"
[1] "Log-like diff: 2.21689333557151e-12"
Running g = 6 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0"
Running g = 7 ...
[1] "Initialization: 1"
[1] "Log-like diff: 1.02318153949454e-12"
Running g = 8 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0"
Running g = 9 ...
[1] "Initialization: 1"
[1] "Log-like diff: 1.01434523003263e-08"
Running g = 10 ...
[1] "Initialization: 1"
[1] "Log-like diff: 2.85737655758567e-10"
Running g = 11 ...
[1] "Initialization: 1"
[1] "Log-like diff: 6.15359752487166e-06"
Running g = 12 ...
[1] "Initialization: 1"
[1] "Log-like diff: 2.24275709115318e-10"
$ICL

$profiles

$boxplots

$probapost_barplots

*************************************************
Model: Poisson
Transformation: none
*************************************************
Clusters fit: 2,3,4,5,6,7,8,9,10,11,12
Clusters with errors: ---
Selected number of clusters via ICL: 12
ICL of selected model: -729818.3
*************************************************
Number of clusters = 12
ICL = -729818.3
*************************************************
Cluster sizes:
 Cluster 1  Cluster 2  Cluster 3  Cluster 4  Cluster 5  Cluster 6  Cluster 7 
        10          9         13         22          3          6          6 
 Cluster 8  Cluster 9 Cluster 10 Cluster 11 Cluster 12 
        19          4         14         23          2 

Number of observations with MAP > 0.90 (% of total):
131 (100%)

Number of observations with MAP > 0.90 per cluster (% of total per cluster):
 Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
 10        9         13        22        3         6         6        
 (100%)    (100%)    (100%)    (100%)    (100%)    (100%)    (100%)   
 Cluster 8 Cluster 9 Cluster 10 Cluster 11 Cluster 12
 19        4         14         23         2         
 (100%)    (100%)    (100%)     (100%)     (100%)    

a <- OntologyProfile(sharedBy3)

a$cluster <- clusteredGenes[a$ensembl_gene_id]

entrezID <- list()
nb_clust = max(clusteredGenes)
for (clust in seq(1:nb_clust)) {
    # print(entrez[entrez$cluster == clust,]$ensembl_transcript_id)
    entrezID[[length(entrezID) + 1]] <- na.omit(a[a$cluster == clust, ]$entrezgene_id)
}
names(entrezID) <- as.character(seq(1:nb_clust))

ck <- compareCluster(geneCluster = entrezID, fun = "enrichGO", OrgDb = org.At.tair.db, 
    ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05)

clusterProfiler::dotplot(ck, x = ~Cluster)

# On essaie un autre clustering avec lka librarie MPLN mpln (dataset =
# as.matrix(data[sharedBy3,])) beaucoup trop long, même mutlithreadé, c'est n'imp

Model-Based Clustering Using MPLN (Parallelized) Description Performs clustering using mixtures of multivariate Poisson-log normal (MPLN) distribution and model selection using AIC, AIC3, BIC and ICL. Since each component/cluster size (G) is independent from another, all Gs in the range to be tested have been parallelized to run on a seperate core using the parallel R package.

Visualisation en ACP

ACP(normalized.count[sharedBy3, ], clusteredGenes)

Class: pca dudi
Call: dudi.pca(df = log(data + 0.1), center = TRUE, scale = TRUE, scannf = FALSE, 
    nf = 4)

Total inertia: 24

Eigenvalues:
    Ax1     Ax2     Ax3     Ax4     Ax5 
17.2861  4.7060  0.6320  0.5028  0.2615 

Projected inertia (%):
    Ax1     Ax2     Ax3     Ax4     Ax5 
 72.025  19.608   2.633   2.095   1.089 

Cumulative projected inertia (%):
    Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
  72.03   91.63   94.27   96.36   97.45 

(Only 5 dimensions (out of 24) are shown)

NULL

Réseau

log.data <- log2(normalized.count[sharedBy3, ] + 1)
Norm.interest.corr <- corr.test(t(log.data), method = "pearson", ci = F)
Norm.interest.corr$p[lower.tri(Norm.interest.corr$p, diag = -TRUE)] = NA
Pval.adj <- as.data.frame(as.table(Norm.interest.corr$p))
Norm.interest.corr$r[lower.tri(Norm.interest.corr$r, diag = TRUE)] = NA
Correlation <- as.data.frame(as.table(Norm.interest.corr$r))
Cor.table <- na.exclude(cbind(Correlation, Pval.adj))[, c(1, 2, 3, 6)]
colnames(Cor.table) <- c("gene1", "gene2", "cor", "p.adj")
Cor.table.filt <- Cor.table[(abs(Cor.table[, 3]) > 0.9 & Cor.table[, 4] < 0.01), 
    ]

g <- graph.data.frame(Cor.table.filt[, 1:2], directed = -FALSE)
V(g)$color <- clusteredGenes[V(g)]
Node_nw_st <- netStats(g)

plot.igraph(g, vertex.size = 5, vertex.label.cex = 0.3, color = clusteredGenes)

# library(d3r) data_json <- d3_igraph(g) write(data_json, 'data.json')

# write.table(Cor.table.filt, 'GraphCO2.txt', sep='\t', row.names=F, quote=F)
# write.table(Node_nw_st,file='StatsCO2.txt', sep='\t', col.names = NA, quote=F)

Réseau avec PLN Network

net_norm <- PLN_network(data = normalized.count, DEGenes = sharedBy3)

 Initialization...
 Adjusting 30 PLN with sparse inverse covariance estimation
    Joint optimization alternating gradient descent and graphical-lasso
    sparsifying penalty = 1.677609 
    sparsifying penalty = 1.549559 
    sparsifying penalty = 1.431282 
    sparsifying penalty = 1.322034 
    sparsifying penalty = 1.221124 
    sparsifying penalty = 1.127917 
    sparsifying penalty = 1.041824 
    sparsifying penalty = 0.9623023 
    sparsifying penalty = 0.8888506 
    sparsifying penalty = 0.8210054 
    sparsifying penalty = 0.7583388 
    sparsifying penalty = 0.7004554 
    sparsifying penalty = 0.6469903 
    sparsifying penalty = 0.597606 
    sparsifying penalty = 0.5519913 
    sparsifying penalty = 0.5098583 
    sparsifying penalty = 0.4709412 
    sparsifying penalty = 0.4349947 
    sparsifying penalty = 0.4017919 
    sparsifying penalty = 0.3711235 
    sparsifying penalty = 0.3427959 
    sparsifying penalty = 0.3166306 
    sparsifying penalty = 0.2924625 
    sparsifying penalty = 0.2701391 
    sparsifying penalty = 0.2495196 
    sparsifying penalty = 0.230474 
    sparsifying penalty = 0.2128821 
    sparsifying penalty = 0.196633 
    sparsifying penalty = 0.1816242 
    sparsifying penalty = 0.1677609 
 Post-treatments
 DONE!

Stability Selection for PLNnetwork: 
subsampling: ++++++++++++++++++++

V(net_norm)$color <- clusteredGenes[V(net_norm)]
plot.igraph(net_norm, vertex.size = 10, vertex.label.cex = 0.5, edge.width = 0.5)

Node_nw_st_pln <- netStats(net_norm)

length(E(net_norm))

[1] 428

length(E(g))

[1] 476

Clustering

load("GenesNitrate_At.RData")
clusteredGenes <- clustering(sharedBy3, round(normalized.count, 0))

****************************************
coseq analysis: Poisson approach & none transformation
K = 2 to 12 
Use set.seed() prior to running coseq for reproducible results.
****************************************
Running g = 2 ...
[1] "Initialization: 1"
[1] "Log-like diff: 6.17148792514399e-07"
Running g = 3 ...
[1] "Initialization: 1"
[1] "Log-like diff: 9.98208899005704e-06"
Running g = 4 ...
[1] "Initialization: 1"
[1] "Log-like diff: 1.6142109870998e-10"
Running g = 5 ...
[1] "Initialization: 1"
[1] "Log-like diff: 3.87143550995006e-07"
Running g = 6 ...
[1] "Initialization: 1"
[1] "Log-like diff: 2.45979485669068e-08"
Running g = 7 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.00150737222311292"
[1] "Log-like diff: 4.44115743505336e-05"
[1] "Log-like diff: 1.4054022159371e-06"
Running g = 8 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.0138274856162859"
[1] "Log-like diff: 0.00232911859524521"
[1] "Log-like diff: 0.000367419357157672"
[1] "Log-like diff: 5.70673573143665e-05"
[1] "Log-like diff: 1.01451924514606e-05"
Running g = 9 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.00455965689279836"
[1] "Log-like diff: 0.00178225206127536"
[1] "Log-like diff: 0.000690376042363994"
[1] "Log-like diff: 0.000266776325485552"
[1] "Log-like diff: 0.000102988618525757"
Running g = 10 ...
[1] "Initialization: 1"
[1] "Log-like diff: 4.08627173626996e-05"
[1] "Log-like diff: 1.21345338399692e-06"
Running g = 11 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.0116560639886529"
[1] "Log-like diff: 0.00308443217002363"
[1] "Log-like diff: 0.000812654937465851"
[1] "Log-like diff: 0.000215029244252207"
[1] "Log-like diff: 5.68554505377961e-05"
Running g = 12 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.743409453403792"
[1] "Log-like diff: 0.725735246714255"
[1] "Log-like diff: 0.647560396066101"
[1] "Log-like diff: 0.607038977642524"
[1] "Log-like diff: 0.503588252926317"
$ICL

$profiles

$boxplots

$probapost_barplots

*************************************************
Model: Poisson
Transformation: none
*************************************************
Clusters fit: 2,3,4,5,6,7,8,9,10,11,12
Clusters with errors: ---
Selected number of clusters via ICL: 12
ICL of selected model: -2967918
*************************************************
Number of clusters = 12
ICL = -2967918
*************************************************
Cluster sizes:
 Cluster 1  Cluster 2  Cluster 3  Cluster 4  Cluster 5  Cluster 6  Cluster 7 
       121        101         27         26        118        157         28 
 Cluster 8  Cluster 9 Cluster 10 Cluster 11 Cluster 12 
        52         35          7         66         99 

Number of observations with MAP > 0.90 (% of total):
830 (99.16%)

Number of observations with MAP > 0.90 per cluster (% of total per cluster):
 Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
 120       101       27        25        118       155       28       
 (99.17%)  (100%)    (100%)    (96.15%)  (100%)    (98.73%)  (100%)   
 Cluster 8 Cluster 9 Cluster 10 Cluster 11 Cluster 12
 51        35        7          64         99        
 (98.08%)  (100%)    (100%)     (96.97%)   (100%)    

a <- OntologyProfile(sharedBy3)

a$cluster <- clusteredGenes[a$ensembl_gene_id]

entrezID <- list()
nb_clust = max(clusteredGenes)
for (clust in seq(1:nb_clust)) {
    # print(entrez[entrez$cluster == clust,]$ensembl_transcript_id)
    entrezID[[length(entrezID) + 1]] <- na.omit(a[a$cluster == clust, ]$entrezgene_id)
}
names(entrezID) <- as.character(seq(1:nb_clust))

ck <- compareCluster(geneCluster = entrezID, fun = "enrichGO", OrgDb = org.At.tair.db, 
    ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05)

clusterProfiler::dotplot(ck, x = ~Cluster)

ACP

ACP(normalized.count[sharedBy3, ], clusteredGenes)

Class: pca dudi
Call: dudi.pca(df = log(data + 0.1), center = TRUE, scale = TRUE, scannf = FALSE, 
    nf = 4)

Total inertia: 24

Eigenvalues:
    Ax1     Ax2     Ax3     Ax4     Ax5 
19.1712  3.3535  0.5249  0.3943  0.1205 

Projected inertia (%):
    Ax1     Ax2     Ax3     Ax4     Ax5 
 79.880  13.973   2.187   1.643   0.502 

Cumulative projected inertia (%):
    Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
  79.88   93.85   96.04   97.68   98.18 

(Only 5 dimensions (out of 24) are shown)

NULL

Relevance network

log.data <- log2(normalized.count[sharedBy3, ] + 1)

Norm.interest.corr <- corr.test(t(log.data), method = "pearson", ci = F)
Norm.interest.corr$p[lower.tri(Norm.interest.corr$p, diag = -TRUE)] = NA
Pval.adj <- as.data.frame(as.table(Norm.interest.corr$p))
Norm.interest.corr$r[lower.tri(Norm.interest.corr$r, diag = TRUE)] = NA
Correlation <- as.data.frame(as.table(Norm.interest.corr$r))
Cor.table <- na.exclude(cbind(Correlation, Pval.adj))[, c(1, 2, 3, 6)]
colnames(Cor.table) <- c("gene1", "gene2", "cor", "p.adj")
Cor.table.filt <- Cor.table[(abs(Cor.table[, 3]) > 0.9 & Cor.table[, 4] < 0.01), 
    ]

g <- graph.data.frame(Cor.table.filt[, 1:2], directed = -FALSE)
V(g)$color <- clusteredGenes[V(g)]

netStats(g)

          degree betweenness Rank_stat
AT4G32950    168    4195.902    684.50
AT2G30660    177    3893.654    682.00
AT5G39130    157    4167.511    680.00
AT1G74590    163    3818.345    679.00
AT5G54370    182    3334.164    677.00
AT3G44300    136    7640.805    673.50
AT3G16180    147    3766.436    671.00
AT1G52070    136    4848.893    669.00
AT5G24070    134    3813.655    659.00
AT2G46680    166    2202.663    658.50
AT4G23400    152    2417.920    658.25
AT2G31880    130    4797.134    657.50
AT1G69490    159    2110.001    654.50
AT2G30670    146    2258.923    653.00
AT1G20180    153    1838.916    644.00
AT2G21880    127    3149.782    643.50
AT5G06730    164    1721.279    643.00
AT5G38900    122    3819.292    639.00
AT1G01680    135    2184.533    638.75
AT5G46890    162    1569.691    633.00
AT1G22530    124    2724.852    632.50
AT5G22270    130    2131.515    629.00
AT4G25250    113    4488.276    628.00
AT4G27400    134    1838.540    627.00
AT1G25400    140    1671.199    626.25
 [ reached 'max' / getOption("max.print") -- omitted 670 rows ]

plot.igraph(g, vertex.size = 5, vertex.label.cex = 0.01, color = clusteredGenes)

PLN Ntework

sample_genes <- sample(sharedBy3, 400)
net <- PLN_network(data = normalized.count, DEGenes = sample_genes)

 Initialization...
 Adjusting 30 PLN with sparse inverse covariance estimation
    Joint optimization alternating gradient descent and graphical-lasso
    sparsifying penalty = 0.9165105 
    sparsifying penalty = 0.8465541 
    sparsifying penalty = 0.7819373 
    sparsifying penalty = 0.7222527 
    sparsifying penalty = 0.6671238 
    sparsifying penalty = 0.6162028 
    sparsifying penalty = 0.5691686 
    sparsifying penalty = 0.5257244 
    sparsifying penalty = 0.4855963 
    sparsifying penalty = 0.4485312 
    sparsifying penalty = 0.4142952 
    sparsifying penalty = 0.3826724 
    sparsifying penalty = 0.3534633 
    sparsifying penalty = 0.3264838 
    sparsifying penalty = 0.3015635 
    sparsifying penalty = 0.2785455 
    sparsifying penalty = 0.2572843 
    sparsifying penalty = 0.237646 
    sparsifying penalty = 0.2195067 
    sparsifying penalty = 0.2027519 
    sparsifying penalty = 0.1872761 
    sparsifying penalty = 0.1729814 
    sparsifying penalty = 0.1597779 
    sparsifying penalty = 0.1475822 
    sparsifying penalty = 0.1363174 
    sparsifying penalty = 0.1259124 
    sparsifying penalty = 0.1163016 
    sparsifying penalty = 0.1074244 
    sparsifying penalty = 0.09922479 
    sparsifying penalty = 0.09165105 
 Post-treatments
 DONE!

Stability Selection for PLNnetwork: 
subsampling: ++++++++++++++++++++

V(net)$color <- clusteredGenes[V(net)]
# V(net)$size <- degree[V(net)]*0.1
plot.igraph(net, vertex.size = 10, vertex.label.cex = 0.4, edge.width = 0.5)

netStats(net)

          degree betweenness Rank_stat
AT5G55110    124   6465.1908    400.00
AT4G19750    102   5682.0143    399.00
AT1G33820     93   5629.1427    398.00
AT4G27400     90   3078.8037    396.50
AT1G65680     73   2966.0060    395.00
AT3G53770     74   2281.3523    394.00
AT4G33710     57   3306.1052    393.50
AT4G12520     58   2114.5319    391.00
AT5G46890     66   1368.4007    390.25
AT5G46900     47   2472.9725    389.50
AT4G11340     63   1247.1850    389.00
AT1G07180     66   1139.8052    388.75
AT5G23990     51   2069.3710    388.00
AT4G36850     40   2415.7477    387.75
AT1G15125     40   1420.4534    385.25
AT1G51840     56    637.2550    384.50
AT1G34047     36   1934.0454    383.50
AT5G24070     53    582.7866    383.00
AT1G51830     55    462.1183    382.50
AT4G15550     41    563.4705    381.00
AT5G20790     37    692.6229    380.75
AT1G10070     33   1161.1155    380.50
AT5G43590     38    360.1270    377.50
AT2G30670     37    347.8894    376.25
AT5G24900     29    326.3185    371.75
 [ reached 'max' / getOption("max.print") -- omitted 375 rows ]

Clustering

load("GenesIron_At.RData")
clusteredGenes <- clustering(sharedBy3, round(normalized.count, 0))

****************************************
coseq analysis: Poisson approach & none transformation
K = 2 to 12 
Use set.seed() prior to running coseq for reproducible results.
****************************************
Running g = 2 ...
[1] "Initialization: 1"
[1] "Log-like diff: 8.36200797493802e-09"
Running g = 3 ...
[1] "Initialization: 1"
[1] "Log-like diff: 4.39287562073507"
[1] "Log-like diff: 2.52480314339526"
[1] "Log-like diff: 0.825138420982821"
[1] "Log-like diff: 0.370149744867518"
[1] "Log-like diff: 0.0851663471061315"
Running g = 4 ...
[1] "Initialization: 1"
[1] "Log-like diff: 88.1138523377031"
[1] "Log-like diff: 84.9471491380293"
[1] "Log-like diff: 154.660538055586"
[1] "Log-like diff: 128.518282504861"
[1] "Log-like diff: 55.9516666470383"
Running g = 5 ...
[1] "Initialization: 1"
[1] "Log-like diff: 148.336500928679"
[1] "Log-like diff: 645.520412267485"
[1] "Log-like diff: 302.729030801015"
[1] "Log-like diff: 241.766725069092"
[1] "Log-like diff: 728.055855867076"
Running g = 6 ...
[1] "Initialization: 1"
[1] "Log-like diff: 3949.59038766225"
[1] "Log-like diff: 107.010732887527"
[1] "Log-like diff: 243.803760768549"
[1] "Log-like diff: 10.2548059812586"
[1] "Log-like diff: 0.617069783586512"
Running g = 7 ...
[1] "Initialization: 1"
[1] "Log-like diff: 923.15734199355"
[1] "Log-like diff: 555.782162059735"
[1] "Log-like diff: 3326.22418578425"
[1] "Log-like diff: 703.479708350589"
[1] "Log-like diff: 1604.5807525359"
Running g = 8 ...
[1] "Initialization: 1"
[1] "Log-like diff: 4.65589725072683"
[1] "Log-like diff: 17.9187198894825"
[1] "Log-like diff: 5.00776285527656"
[1] "Log-like diff: 0.37400244270319"
[1] "Log-like diff: 0.049682970164266"
Running g = 9 ...
[1] "Initialization: 1"
[1] "Log-like diff: 416.347468761451"
[1] "Log-like diff: 323.722501753506"
[1] "Log-like diff: 499.810454781299"
[1] "Log-like diff: 932.723909985434"
[1] "Log-like diff: 414.498517703427"
Running g = 10 ...
[1] "Initialization: 1"
[1] "Log-like diff: 508.448868045691"
[1] "Log-like diff: 174.336283681083"
[1] "Log-like diff: 92.140530301998"
[1] "Log-like diff: 764.239944442989"
[1] "Log-like diff: 1952.41368072682"
Running g = 11 ...
[1] "Initialization: 1"
[1] "Log-like diff: 749.173439262413"
[1] "Log-like diff: 659.603164162899"
[1] "Log-like diff: 1673.26394014634"
[1] "Log-like diff: 672.267399303077"
[1] "Log-like diff: 770.792673019672"
Running g = 12 ...
[1] "Initialization: 1"
[1] "Log-like diff: 59.9308044587635"
[1] "Log-like diff: 33.9311456470273"
[1] "Log-like diff: 17.1399382118373"
[1] "Log-like diff: 5.13560451611352"
[1] "Log-like diff: 2.11074602623562"
$ICL

$profiles

$boxplots

$probapost_barplots

*************************************************
Model: Poisson
Transformation: none
*************************************************
Clusters fit: 2,3,4,5,6,7,8,9,10,11,12
Clusters with errors: ---
Selected number of clusters via ICL: 11
ICL of selected model: -3313489
*************************************************
Number of clusters = 11
ICL = -3313489
*************************************************
Cluster sizes:
 Cluster 1  Cluster 2  Cluster 3  Cluster 4  Cluster 5  Cluster 6  Cluster 7 
       714        218         31        106        411        594        221 
 Cluster 8  Cluster 9 Cluster 10 Cluster 11 
        15        294         76        161 

Number of observations with MAP > 0.90 (% of total):
2759 (97.11%)

Number of observations with MAP > 0.90 per cluster (% of total per cluster):
 Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
 704       204       29        101       398       588       209      
 (98.6%)   (93.58%)  (93.55%)  (95.28%)  (96.84%)  (98.99%)  (94.57%) 
 Cluster 8 Cluster 9 Cluster 10 Cluster 11
 13        284       73         156       
 (86.67%)  (96.6%)   (96.05%)   (96.89%)  

a <- OntologyProfile(sharedBy3)

a$cluster <- clusteredGenes[a$ensembl_gene_id]

entrezID <- list()
nb_clust = max(clusteredGenes)
for (clust in seq(1:nb_clust)) {
    # print(entrez[entrez$cluster == clust,]$ensembl_transcript_id)
    entrezID[[length(entrezID) + 1]] <- na.omit(a[a$cluster == clust, ]$entrezgene_id)
}
names(entrezID) <- as.character(seq(1:nb_clust))

ck <- compareCluster(geneCluster = entrezID, fun = "enrichGO", OrgDb = org.At.tair.db, 
    ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05)

clusterProfiler::dotplot(ck, x = ~Cluster)

ACP

ACP(normalized.count[sharedBy3, ], clusteredGenes)

Class: pca dudi
Call: dudi.pca(df = log(data + 0.1), center = TRUE, scale = TRUE, scannf = FALSE, 
    nf = 4)

Total inertia: 24

Eigenvalues:
     Ax1      Ax2      Ax3      Ax4      Ax5 
22.01676  1.12356  0.26987  0.11457  0.06979 

Projected inertia (%):
    Ax1     Ax2     Ax3     Ax4     Ax5 
91.7365  4.6815  1.1245  0.4774  0.2908 

Cumulative projected inertia (%):
    Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
  91.74   96.42   97.54   98.02   98.31 

(Only 5 dimensions (out of 24) are shown)

NULL

Relevance network

sample_genes <- sample(sharedBy3, 400)
log.data <- log2(normalized.count[sharedBy3, ] + 1)


Norm.interest.corr <- corr.test(t(log.data), method = "pearson", ci = F)
Norm.interest.corr$p[lower.tri(Norm.interest.corr$p, diag = -TRUE)] = NA
Pval.adj <- as.data.frame(as.table(Norm.interest.corr$p))
Norm.interest.corr$r[lower.tri(Norm.interest.corr$r, diag = TRUE)] = NA
Correlation <- as.data.frame(as.table(Norm.interest.corr$r))
Cor.table <- na.exclude(cbind(Correlation, Pval.adj))[, c(1, 2, 3, 6)]
colnames(Cor.table) <- c("gene1", "gene2", "cor", "p.adj")
Cor.table.filt <- Cor.table[(abs(Cor.table[, 3]) > 0.9 & Cor.table[, 4] < 0.01), 
    ]

g <- graph.data.frame(Cor.table.filt[, 1:2], directed = -FALSE)
V(g)$color <- clusteredGenes[V(g)]


netStats(g)

          degree betweenness Rank_stat
AT3G24750    408    29895.80   2422.50
AT3G59660    374    30464.74   2421.50
AT4G16920    464    24536.84   2418.00
AT5G55135    333    24067.60   2402.25
AT4G04490    385    16625.79   2392.50
AT1G49160    312    22683.14   2391.00
AT4G15563    304    23216.59   2388.75
AT1G77330    335    17580.56   2384.75
AT5G56320    349    14437.99   2373.75
AT3G62120    253    50151.01   2369.25
AT2G07050    246    64836.94   2368.50
AT1G09935    272    21354.01   2367.00
AT3G17770    332    13362.11   2359.50
AT2G17520    244    32437.83   2359.00
AT1G69030    291    15125.51   2358.25
AT5G51640    322    13584.21   2357.00
AT2G17220    371    12281.28   2356.75
AT5G51070    364    12367.83   2355.50
AT3G17700    370    12199.02   2354.00
AT5G24070    375    11794.52   2349.50
AT5G65810    297    13433.66   2347.75
AT5G19760    255    19783.74   2347.75
AT5G53350    223    40735.04   2346.25
AT2G31880    257    18948.39   2345.75
AT5G12900    280    14435.61   2345.25
 [ reached 'max' / getOption("max.print") -- omitted 2411 rows ]

plot.igraph(g, vertex.size = 5, vertex.label.cex = 0.01, color = clusteredGenes)

PLN Network

PLN_network(data = data, DEGenes = sample_genes)

 Initialization...
 Adjusting 30 PLN with sparse inverse covariance estimation
    Joint optimization alternating gradient descent and graphical-lasso
    sparsifying penalty = 1.646713 
    sparsifying penalty = 1.521021 
    sparsifying penalty = 1.404923 
    sparsifying penalty = 1.297686 
    sparsifying penalty = 1.198635 
    sparsifying penalty = 1.107144 
    sparsifying penalty = 1.022637 
    sparsifying penalty = 0.9445798 
    sparsifying penalty = 0.8724808 
    sparsifying penalty = 0.8058851 
    sparsifying penalty = 0.7443726 
    sparsifying penalty = 0.6875553 
    sparsifying penalty = 0.6350748 
    sparsifying penalty = 0.5866001 
    sparsifying penalty = 0.5418254 
    sparsifying penalty = 0.5004683 
    sparsifying penalty = 0.462268 
    sparsifying penalty = 0.4269835 
    sparsifying penalty = 0.3943922 
    sparsifying penalty = 0.3642886 
    sparsifying penalty = 0.3364827 
    sparsifying penalty = 0.3107993 
    sparsifying penalty = 0.2870763 
    sparsifying penalty = 0.265164 
    sparsifying penalty = 0.2449242 
    sparsifying penalty = 0.2262294 
    sparsifying penalty = 0.2089615 
    sparsifying penalty = 0.1930116 
    sparsifying penalty = 0.1782792 
    sparsifying penalty = 0.1646713 
 Post-treatments
 DONE!

Stability Selection for PLNnetwork: 
subsampling: ++++++++++++++++++++

IGRAPH e8659a4 UNW- 400 365 -- 
+ attr: name (v/c), label (v/c), label.cex (v/n), size (v/n),
| label.color (v/c), frame.color (v/l), weight (e/n), color (e/c),
| width (e/n)
+ edges from e8659a4 (vertex names):
 [1] AT4G24040--AT4G33560 AT1G15640--AT3G22235 AT1G15640--AT1G08100
 [4] AT1G15640--AT1G33820 AT4G30640--AT1G08100 AT1G75890--AT3G22235
 [7] AT1G75890--AT2G36970 AT1G75890--AT1G08630 AT5G35940--AT1G08100
[10] AT4G33560--AT3G43850 AT4G33560--AT5G24550 AT4G33560--AT5G25840
[13] AT4G33560--AT1G31772 AT4G33560--AT1G31580 AT4G33560--AT1G62500
[16] AT4G33560--AT1G20990 AT4G33560--AT5G15120 AT5G02920--AT3G46230
+ ... omitted several edges

net <- PLN_network(data = normalized.count, DEGenes = sample_genes)

 Initialization...
 Adjusting 30 PLN with sparse inverse covariance estimation
    Joint optimization alternating gradient descent and graphical-lasso
    sparsifying penalty = 1.64184 
    sparsifying penalty = 1.51652 
    sparsifying penalty = 1.400765 
    sparsifying penalty = 1.293846 
    sparsifying penalty = 1.195088 
    sparsifying penalty = 1.103868 
    sparsifying penalty = 1.01961 
    sparsifying penalty = 0.9417844 
    sparsifying penalty = 0.8698988 
    sparsifying penalty = 0.8035002 
    sparsifying penalty = 0.7421697 
    sparsifying penalty = 0.6855205 
    sparsifying penalty = 0.6331954 
    sparsifying penalty = 0.5848641 
    sparsifying penalty = 0.5402219 
    sparsifying penalty = 0.4989872 
    sparsifying penalty = 0.4609 
    sparsifying penalty = 0.4257199 
    sparsifying penalty = 0.393225 
    sparsifying penalty = 0.3632105 
    sparsifying penalty = 0.335487 
    sparsifying penalty = 0.3098795 
    sparsifying penalty = 0.2862267 
    sparsifying penalty = 0.2643793 
    sparsifying penalty = 0.2441994 
    sparsifying penalty = 0.2255599 
    sparsifying penalty = 0.2083431 
    sparsifying penalty = 0.1924404 
    sparsifying penalty = 0.1777516 
    sparsifying penalty = 0.164184 
 Post-treatments
 DONE!

Stability Selection for PLNnetwork: 
subsampling: ++++++++++++++++++++

V(net)$color <- clusteredGenes[V(net)]
# V(net)$size <- degree[V(net)]*0.1

plot.igraph(net, vertex.size = 10, vertex.label.cex = 0.4, edge.width = 0.5)

netStats(net)

          degree betweenness Rank_stat
AT1G08100    120  7672.01950    400.00
AT1G33820     55  1454.87124    399.00
AT3G22235     36   857.01867    398.00
AT2G29350     33   324.99423    397.00
AT5G25840     22   201.49699    395.00
AT1G08630     18   266.09726    394.75
AT5G24550     12   158.81180    391.75
AT1G54970     18    92.23680    391.75
AT4G23140     13   129.60329    391.50
AT4G33560      9   293.57083    390.25
AT1G31580     11   134.75842    390.25
AT4G15550     12    88.86616    389.25
AT5G10760     15    58.34160    388.50
AT1G31772      9    96.72239    387.25
AT2G36970     11    63.74967    387.25
AT3G43850      9    62.67112    384.75
AT5G46610      9    36.60439    381.75
AT3G13433     10    22.12647    381.00
AT3G46230      6    43.97498    378.50
AT1G66600      8    23.54934    378.50
AT5G35810      7    25.70546    378.00
AT1G63550      8    21.00140    377.50
AT4G31950      6    42.51471    377.50
AT1G62500      5    71.42768    377.00
AT4G36820      6    28.47573    376.50
 [ reached 'max' / getOption("max.print") -- omitted 375 rows ]

Enquête sur la similarité entre cNF et CNF

• Corrélations entre les réplicats à l’intérieur d’une condition et de l’autre sont faiblement supérieures à celles entre cNF et CNF (images investigations factorCO2)

Réseau CO2

• Quand on compare CNF à une condition x et cNF à cette même condition x (6 conditions possibles pour x), on retrouve entre 40 et 70% de gènes en commun, suggérant toute fois des différences entre ces deux transcritômes (on aurait presque 100% de similarité sinon)

Application des ces méthodes à la tomate

• La tomate semble répondre différemment dans certaines mesures : plus d’effet du CO2, effet moindre du fer, effet nitrate plutôt similaire.

• Ontologies moins fournies pour la tomate

Réseau CO2

Relevance network sur les gènes qui répondent globalement à un facteur

• DEG en commun entre les 4 comparaisons possibles pour l’effet d’un facteur (Venn diagrams)

• Fait sur l’ensemble des transcriptômes (plus large que les transcriptômes sur lesquels les DEG ont été détectés)

• Relevance Network fait comme Rodrigo. et Al, seuil sur la valeur de corréation et sur la pvalue, gènes triés sur leur centralité et connectivité

• Visualisés dans igraph après Clustering

• Visualisés dans Cytoscape, pus clustering de communautés (pluggins) et analyse d’enrichissement d’ontologies

Réseau CO2

Début de biblio sur les méthodes d’inférence de GRN

• Regression: Transcription factors are selected by target gene specific sparse linear regression and data resampling approaches.

• Bayesian networks optimize posterior probabilities by different heuristic searches.

• Correlation Edges are ranked based on variants of correlation.

• Méthodes mixtes Consensus entre différentes techniques (conclusion de la métanalyse et benchmarcks du projet DREAM5 (wisdom of crowds))

Others : Genie3: A random forest, neural networks, chi 2 - Mutual Information: Edges are (1) ranked based on variants of mutual information and (2) filtered for causal relationships.